﻿using System;
using System.Collections.Generic;
using System.Diagnostics.CodeAnalysis;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

using MathNet.Numerics.LinearAlgebra;
using MathNet.Numerics.LinearAlgebra.Double;
using MathNet.Numerics.Statistics;
using MathNet.Symbolics;

using Expr = MathNet.Symbolics.SymbolicExpression;


namespace NumericalAnalysis.DataInterpolation
{
    public class AitkenInterpolation
    {
        Expr polynomial = 0;   //多项式表达式
        List<double> polyCoefficients = new List<double>(); //各项系数 幂次高至低
        List<int> coefficienOrders = new List<int>(); //各项幂次

        Vector<double> theVec = Vector<double>.Build.Dense(2);
        List<Vector<double>> vectors = new List<Vector<double>>();

        int count = 0;

        public AitkenInterpolation(IEnumerable<Vector<double>> _vectors)
        {
            count = _vectors.Count();
            vectors = _vectors.ToList();
        }
        public AitkenInterpolation(IEnumerable<Tuple<double, double>> _tuples)
        {
            count = _tuples.Count();
            vectors = _tuples.Select(v =>
            {
                var vec = Vector<double>.Build.Dense(2);
                vec[0] = v.Item1;
                vec[1] = v.Item2;
                return vec;
            }).ToList();
        }

        /// <summary>
        /// Accuracy most
        /// </summary>
        /// <returns></returns>
        public Expr FitInterpolation()
        { 
            if (polynomial.ToString().Contains("t")) return polynomial;
            var t = Expr.Variable("t");

            List<Expr> polyNexts = vectors.Select(v => (Expr)0).ToList();
            List<Expr> polyPervs = vectors.Select(v => (Expr)v[1]).ToList();
            //TODO  迭代算法
            for (int i = 0; i < count - 1; i++)
            {
                for (int j = i + 1; j < count; j++)
                {
                    var xk = vectors[i][0];
                    var xn = vectors[j][0];
                    var pk = polyPervs[i];
                    var pn = polyPervs[j];
                    //var matrix =  Matrix<Expr>.Build.DenseIdentity(2);
                    polyNexts[j] = (pk * (xn - t) - pn * (xk - t)) / (xn - xk);
                    polyNexts[j] = polyNexts[j].Expand();
                }
                polyPervs = polyNexts;
            }
            polynomial = polyNexts.Last();


            //显示多项式
            Console.WriteLine("--------------------------------");
            Console.WriteLine("多项式 polynomial:");
            Console.WriteLine(polynomial);

            //Console.WriteLine("--------------------------------");
            //Console.WriteLine("多项式简化 polynomial-ExponentialSimplify:");
            //Console.WriteLine(polynomial.ExponentialSimplify());

            //Console.WriteLine("--------------------------------");
            //Console.WriteLine("多项式展开式 polynomial-Expand:");
            //Console.WriteLine(polynomial.Expand());

            polynomial = polynomial.Expand();

            return polynomial;
        }

        /// <summary>
        /// 从给出的数列计算
        /// </summary>
        /// <param name="xValues"></param>
        /// <returns></returns>
        public IList<double> CalculateBy(IList<double> xValues)
        {
            List<double> yValues = new List<double>();
            for (int i = 0; i < xValues.Count(); i++)
            {
                var symbols = new Dictionary<string, FloatingPoint> { { "t", xValues[i] } };
                double yValue = polynomial.Evaluate(symbols).RealValue;

                yValues.Add(yValue);
            }
            return yValues;
        }

        /// <summary>
        /// 从给出的数计算
        /// </summary>
        /// <param name="xValue"></param>
        /// <returns></returns>
        public double CalculateBy(double xValue)
        {
            //复杂多项式避免使用Compile方法计算  将严重拖慢性能
            //Func<double, double> func = polynomial.Compile("t"); 
            //double yValue = func(xValue);

            var symbols = new Dictionary<string, FloatingPoint> { { "t", xValue } };
            double yValue = polynomial.Evaluate(symbols).RealValue;

            return yValue;
        }

    }


    public class AitkenInterpolationTest
    {
        /// <summary>
        /// 测试
        /// </summary>
        /// <param name="exsistPnts">存在插值点</param>
        /// <param name="needPntXs">需要求的点自变量</param>
        /// <returns></returns>
        public static IList<Tuple<double, double>> AitkenInterpolationCalc(IList<Tuple<double, double>> exsistPnts, IList<double> needPntXs)
        {
            var t = Expr.Variable("t");

            AitkenInterpolation lagInterp = new AitkenInterpolation(exsistPnts);
            var expression = lagInterp.FitInterpolation();

            var needPntYs = lagInterp.CalculateBy(needPntXs);

            var tuples = new List<Tuple<double, double>>();
            for (int i = 0; i < needPntXs.Count; i++)
            {
                var tuple = new Tuple<double, double>(needPntXs[i], needPntYs[i]);
                tuples.Add(tuple);
            }
            return tuples;
        }



        /// <summary>
        /// sin函数
        /// </summary>
        /// <returns></returns>
        public static (List<Tuple<double, double>>, List<Tuple<double, double>>) ExampleSin()
        {
            var t = Expr.Variable("t");
            List<Tuple<double, double>> extVecs = new List<Tuple<double, double>>();
            //构造sin函数插值点
            List<double> array = NumericalService.LinSpace(0, Math.PI * 2, Math.PI * 2 * 0.1);
            extVecs = NumericalService.CalcByFunc((x) => { return Math.Sin(x); }, array).ToList();

            //进行插值 并 返回检查点    这里插值点就是检查点
            var interpYs = AitkenInterpolationCalc(extVecs, extVecs.Select(v => v.Item1).ToList()).ToList();

            List<double> array2 = NumericalService.LinSpace(0 + 0.5, 9 - 0.5, 0.01);
            var checkdYs = AitkenInterpolationCalc(extVecs, array2).ToList();

            return (interpYs, checkdYs);

        }



    }

}
